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Abstract 

This paper presents a counterexample to the conjecture that the semi-expHcit Lie- 
Newmark algorithm is variational. As a consequence the Lie-Newmark method is not 
well-suited for long-time simulation of rigid body-type mechanical systems. The coun- 
terexample consists of a single rigid body in a static potential field, and can serve as a test 
of the variational nature of other rigid-body integrators. 
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1 Introduction 

In this paper we will focus on the dynamics of a rigid body in a static potential field. To 
describe this system, denote by Q{t) G S0(3), W{t) = [Wi{t) W2{t) W3{t)]^ G R\ and I = 
diag(/i,/2,/3) G M^'^^ the configuration, body angular velocity and inertia matrix of the body, 
respectively. Let T : S0(3) — )■ M? be the torque acting on the body and ^: M? — )• M.^^^ be the 
hat map 

-W3 W2 ' 
W= IV3 -Wi 
-W2 Wi 



In terms of this notation, the governing equations are 

[q =qw 

\ m =IW xW + x{Q), 



(la) 
(lb) 



with initial conditions g(0) = 2o G S0(3) and W{Q) = Wq G M^. We assume this rigid body is 
derivable from a Lagrangian L : S0(3) x M^ — )• M of the form 



L{Q,W) = T(W)-U{Q), 



(2) 
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where T{W) = ^W'^IW and U{Q) are the kinetic and potential energy of the body, respectively. 
This assumption implies that the torque in (lb) can be computed from the directional derivative 
of L'^ at 2 in the direction Qy: 

ziQfy = -DU{Q)-Qy, yGR\ 

Notice that the total energy is separable and T{—W) = T{W). The exact flow of (1) possesses 
certain structure such as total energy preservation, time-symmetry, and symplecticity. More- 
over, the path Q lies on a configuration manifold SO (3) which possesses a Lie-group structure. 

This paper investigates the long-run behavior of two integrators for (1): the 
Lie-Newmark [SVQ88, SW91] and Lie-Verlet methods [BRM08]. Both methods are semi- 
explicit, second-order accurate and symmetric. They are also 'Lie group methods' because 
they respect the Lie group structure of the configuration manifold [IMKNZOO]. Moreover, the 
complexity and implementation of the two methods is quite similar. The main difference be- 
tween the integrators is that the Lie-Verlet method is designed to be variational, whereas the 
Lie-Newmark method is not. 

Variational integrators are time-integrators adapted to the structure of mechanical systems 
[MWOl, LMOW04a, LMOW04b]. They are symplectic, and in the presence of symmetry, 
momentum preserving. The theory of variational integrators includes discrete analogues of 
Hamilton's principle, Noether's theorem, the Euler-Lagrange equations, and the Legendre trans- 
form. The variational nature of Lie-Verlet guarantees its excellent long-time behavior. In 
fact, one can prove this. The basic idea of the proof is to show that a trajectory of a varia- 
tional integrator is interpolated by a level set of a 'modified' energy function nearby the true 
energy [BG94, Rei99, HLW06]. This implies that a trajectory of the variational integrator is 
confined to these level sets for the duration of the simulation. As a consequence variational 
integrators nearly preserve the true energy and exhibit linear growth in global error. For these 
reasons variational integrators are well-suited for long-time simulation. 

Even though the Lie-Newmark integrator is not designed to be variational, this does not 
rule out the possibility that the algorithm is variational in a subtle way like classical Newmark. 
The classical Newmark family of algorithms are widely used integrators in computational me- 
chanics (for an expository treatment see, e.g.. Chapter 9 of [HugSV]). They were first proposed 
in [New59], but it took more than forty years for their variational nature to be established 
in [KMOWOO]. Specifically, Kane et al. proved that a trajectory of the classical Newmark 
method is shadowed by a trajectory of a variational algorithm. In other words the Newmark 
integrators are not symplectic, but conjugate symplectic [HLW06]. The possibility that Lie- 
Newmark could be analogously variational (or conjugate symplectic) was supported by recent 
numerical evidence showing that the Lie-Newmark algorithm exhibits excellent behavior akin 
to classical Newmark [KE05]. Based on these experiments, Krysl and Endres conjectured that 
the Lie-Newmark algorithm is variational. 

This paper disproves this conjecture. In particular, the paper presents a simple numerical 
counterexample showing that the Lie-Newmark method exhibits systematic energy drift. In 
contrast, the Lie-Verlet method nearly preserves the true energy and exhibits the qualitative 
properties one expects of a variational integrator. In summary, the Lie-Verlet method is well- 
suited for long-time simulation of rigid body-type mechanical systems, while the Lie-Newmark 
method is not. 



The paper is organized as follows. The algorithms used in the paper are stated in §2 and the 
numerical 'stress test' carried out in §3. In §4 we provide concluding remarks including other 
applications of this numerical 'stress test.' Along with this paper, we have released a simple 
Matlab implementation of the presented algorithms. This release can be retrieved from the 
Matlab Central File Exchange. 

2 Integrators 

Lie-Newmark The Lie-Newmark family of integrators was proposed more than 
twenty years ago in [SVQ88]. These methods consist of a Newmark-style discretization of (lb) 
and a discretization of (la) that ensures the configuration update remains on SO (3). These meth- 
ods were motivated by the need to develop conserving algorithms that can efficiently simulate 
the structural dynamics of rods and shells. For example, consider simulating large deformations 
of a three-dimensional finite-strain rod model. The rod is typically discretized using N copies 
of M.^ X S0(3) where N is the number of discretization points along the line of centroids of the 
rod. The configuration of the rod at each point along the line of centroids is specified by an 
element of S0(3). The dynamical behavior of the rod can then be estimated by simulating the 
dynamics of A'^ rigid bodies with torques due to elastic coupling between bodies. 

This paper focuses on a specific member of the Lie-Newmark family tested in [KE05]. This 
member is the Lie group analog of the so-called explicit Newmark method on vector spaces 
(see Chapter 9 of [Hug87]). In molecular dynamics the explicit Newmark method is known 
as the Verlet integrator [HLW06]. Given {Qk,Wk) G M^ x S0(3) and time-stepsize h, the Lie- 
Newmark algorithm determines (2yt+i, Wyt+i) by the following iteration rule: 

Wi^^^^=Wk + '^r\lWkXWk + TiQk)), (3a) 

Qk+i=Qkcay{hW^^l^), (3b) 

Wk+i = W,^i + h-\lWk+, X Wk+i + T(e,+i)) . (3c) 

Here we have introduced the Cay ley map cay : R-' — )■ S0(3): 

4 2 

cay (x) = / H —x ^ —o?- , x G M^ , (4) 

4 + |x| 4 + |x| 

where / is the identity matrix. The Cayley map is a second-order approximation of the expo- 
nential map on S0(3). There are other maps one can use in place of the Cayley map in (3b) 
(see, e.g., [BR07, §5.4]), but the Cayley map is known to be very computationally efficient in 
practice. This integrator is semi-explicit because (3a)-(3b) involve explicit updates, and (3c) 
is only implicit in the angular velocity and not in the torque. Hence, the implicitness of the 
Lie-Newmark method is not severe. It is also symmetric and second-order accurate. 

Lie- Verlet The Lie- Verlet integrator was proposed in [BRM08] and based on the theory of 
discrete and continuous Euler-Poincare systems [MPS98,MS93]. The method is closely related 
to, but different from the RATTLE method for constrained mechanical systems [HLW06]. 



Given {Qk,Wk) € M? x S0(3) and time-stepsize h, the Lie-Verlet algorithm determines 
{Qk+i,Wic+i) by the following iteration rule: 



W, 



k+{ 



-wu + h-' 



IW.^ixW, 



h 



'k+ 



k+i 



C'lw^.+i w^.+i + <a 



Qk+1 =Qkcay{hW,^^ 



Wk+i=w,+i+'^r' iw, 



^+ixW^+i + 



2V ^^ 



.iw,+.)w,+.+T(e,+i, 



(5a) 
(5b) 

(5c) 



Similar to the Lie-Newmark method, this algorithm is symmetric, semi-explicit and second- 
order accurate. In particular, the updates in (5b) and (5c) are explicit, and the implicitness in 
(5a) does not involve the torque. We emphasize the Lie-Verlet integrator is variational and refer 
the reader to [BRM08] for a proof of this result. 

3 Numerical Counterexample 

This section describes a numerical experiment showing that the semi-explicit Lie-Newmark in- 
tegrator (3) exhibits systematic drift in total energy. Such drift implies that the method is not a 
conjugate symplectic integrator for (1), and hence, is not variational. The numerical counterex- 
ample we discuss is strongly inspired by a numerical experiment reported in [FHP04, §4.4]. 
That paper shows systematic energy drift along an orbit of a fourth-order accurate, implicit, and 
symmetric Lobatto IIIB integrator when applied to a spring pendulum with exterior forces. 



Preliminaries Let / G S0(3) be the identity matrix, and in terms of which, define the function 
dist : S0(3) x S0(3) ^ M as 



dist(ei,e2 



2tr(/-e[e2). 



(6) 



Let II • \\f denote the Frobenius matrix norm. We recall that \\A\\f := A/tr(A^A) for A G M"^". 
It is straightforward to verify that dist(-, •) is a metric on S0(3) induced by the Frobenius norm 
using the identity 

IIG2-eiiiF = 2tr(/-e[e2). 

For the numerical experiment, consider a single rigid body in a static field defined by the 
potential energy function Ua '■ S0(3) — )■ M given by: 

a 



Ua{Q) = {dht{Q,I)-\) 



dist(e,e» 



(7) 



The first term in the right hand side of (7) is a bounded potential which attains its minimum 
value at 2 G S0(3) satisfying dist(2,/) = 1. The second term is an unbounded potential that 
generates an attraction toward the configuration Q,„ G S0(3). The parameter a is a tuning 
parameter. 

For a = 0, the potential energy Ua=o achieves its minimum value on the two-dimensional 
surface 

5 = {eGSO(3):dist(e,/) = l}. 



This implies that the set S x {0} C S0(3) x M^ is a (locally) stable set in the sense of Lyapunov 
for the dynamics of the rigid body. One proves this fact using as Lyapunov function the energy 

EiQ,W) = T{W) + Ua=o{Q) 

and noting that the set {iQ,W) £ S0(3) x R^ \ E{Q,W) < E} is compact for every E>0. 

For a > 0, the set S gets perturbed by the unbounded attractive potential. On this perturbed 
energy landscape, the rigid body experiences an attraction toward the configuration Q,„. Yet, if 
we place the attraction point Q^ sufficiently far from the set S and choose the tuning parameter 
a > sufficiently small, the set S gets only slightly perturbed into a new set, that we label Sa- 
Furthermore, the set Sa x {0} C S0(3) x M^ is locally Lyapunov stable like the unperturbed set 
Sx{0}. 

In summary, we can design Ua so that the true solution conserves a compact energy function 
when initialized in a neighborhood of the set Sa x {0}. This conserved quantity implies that the 
true solution is confined to this neighborhood. This property will be preserved by the Lie-Verlet 
integrator since it is variational. Recall that a variational integrator is interpolated by a level set 
of a 'modified' energy function nearby the true energy [BG94, Rei99,HLW06]. This implies 
that a trajectory of the variational integrator is confined to a neighborhood of Sa x {0} for the 
duration of the simulation. Next we show that the Lie-Newmark integrator exhibits systematic 
drift in this energy function. 



Numerical 'Stress Test' Now we are in position to describe the numerical counterexample. 
Select the inertia matrix to be 

I: 

and the tuning parameter to be a = 0.3. Place the attraction point at 

2m=exp(Vm), Vm 



'2 





01 





2 











4 



1 

V2 



2.5 



2.5 



Here exp : 



S0(3) is the exponential map 



exp (x) = / + 



sm X 



-x+- 



1 — cosfUl 



xG 



The initial condition (2o,W()) € M^ x S0(3) is selected so that dist(2o5^) is nearly one. Specif- 
ically, select the initial configuration and angular velocity to be 



2o = exp(vo), vo 



In the numerical experiment we test the two integrators, Lie-Newmark (NMB) and velocity 
Lie-Verlet (VLV), on a long time interval [0, 15000]. The energy error obtained with a time- 
stepsize h = 0.125 is shown in Figure 1(a). The experiment was repeated with a time-stepsize 











0.7227 


, Wo = 










0.625 



h = 0.25 and results are reported in Figure 1(b). A systematic drift for the NMB scheme can 
be observed in both cases. The drift appears linear in the time span T and quadratic in the time- 
stepsize h. We abbreviate this fact by saying the total energy error behaves like &{Th^). No 
energy drift is observed for the VLV scheme. We have also tested the Lie-Newmark method 
with the Cayley map in (3b) replaced by the exponential map. Systematic energy drift is ob- 
served in that case too. 
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Figure 1 : This figure shows the energy error of the Lie-Newmark (NMB) and velocity Lie- 
Verlet (VLV) algorithms for the rigid body in the potential energy landscape defined by (7) 
for two different timesteps. NMB exhibits a systematic energy drift. On the other hand, the 
energy error of VLV method remains bounded as predicted by theory. The initial conditions 
and parameters used are provided in the text. 



The trajectory generated by Lie-Newmark for time-stepsize h = 0.25 is shown in the 
axis/angle representation of S0(3) in Figure 2. The semi-transparent surfaces correspond to 
isosurfaces of the potential energy (7). 

The time-precision diagrams, shown in Figures 3(a) and 3(b) confirm that NMB and VLV 
are second-order accurate. Observe from the figures that the slope of the two lines denoting the 
global error is ^(/j^). The diagrams have been generated by computing the global error in the 
configuration and angular velocity evaluated at T = 5. The simulations have been performed 
for a variety of time-stepsizes as indicated in the figures. The reference solution was computed 
using the function ode45 in Matlab, with an absolute tolerance 10^ ''^ and relative tolerance 
2-10-^1 



4 Conclusion 



The Lie-Newmark method was proposed as a generalization of the explicit Newmark algorithm 
to Lie groups [SVQ88]. However, unlike its counterpart on vector spaces, this paper shows 
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Figure 2: This figure shows a trajectory generated by the Lie-Newmark integrator operated at 
time-stepsize h = 0.25 projected on the axis/angle representation of S0(3). The initial condi- 
tions and parameters used are provided in the text. The semi-transparent surfaces are level sets 
of the potential energy (7). The darker shading corresponds to lower potential energy. The dot 
in the figure corresponds to the attraction point Q^ G S0(3) of the potential energy. 

that the Lie-Newmark method does not possess excellent long-time behavior when applied to 
a rigid body in a potential force field. In particular, the paper presents a numerical experiment 
showing the Lie-Newmark integrator exhibits systematic energy drift with an ^'{Th^) behav- 
ior. In contrast, the Lie-Verlet method, which is variational by construction, does not exhibit 
energy drift as predicted by theory. Since the two methods are semi-explicit, symmetric and 
computationally similar to implement, we conclude that the Lie-Verlet method is better suited 
for long-time simulation of rigid body-type systems. 

As a final remark, we note that the numerical example appearing in this paper can serve as 
a test for symplecticity of other rigid body integrators. For instance, in Appendix A we show 
non-symplecticity of the 'explicit' Lie-Midpoint algorithm proposed in [Kry05]. 
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Figure 3: This figure shows the global error of the Lie-Newmark (NMB) and velocity Lie- 
Verlet (VLV) algorithms. The global error is evaluated in configuration and body angular 
velocity at a physical time of T = 5 for a variety of time-stepsizes. We use as a reference 
solution an integration of (1) using the Matlab function ode45 with low tolerance. Observe 
that both integrators are second-order accurate. 



Appendix A Non-Symplecticity of Explicit Lie-Midpoint Algorithm 

In this section the LIEMID[EA] algorithm is subjected to the numerical 'stress test' described in 
§3. This algorithm was thought to be a symplectic-momentum integrator in [Kry05]. However, 
the test reveals that the integrator is not conjugate symplectic. 

Given {Qk.'^k) G K^ x S0(3) and time-stepsize h, the LIEMID[EA] algorithm determines 
(QA+i)W^/t+i) by the following iteration rule: 







k+\ 



Viexp(-i0,+ ,)(W, + ^T(a)), 



G,+ .=e-texp(0, 



k+l,. 



w, 



k+i 



--\ 



-1 



exp(-0,+ .)(W, + -T(e;t)), 



^riexp(-^0,+i)(w,+ .; 



0^+1 

G/t+i = 2«:+iexp(0i+i 



W,+i =r^(exp(-0,)IW,+ . + -T(G^+i)) . 



(8a) 
(8b) 

(8c) 

(8d) 
(8e) 

(8f) 



This integrator is more implicit than the semi-explicit Lie-Newmark and Lie-Verlet methods 
introduced in §2. In particular, the updates (8a) and (8d) are both implicit. Hence, the algorithm 



involves two nonlinear solves per step, unlike the Lie-Newmark and Lie-Verlet methods which 
involve one nonlinear solve per step. The four remaining updates are explicit. This algorithm 
was derived as a composition of a half-step of a first-order Lie-midpoint method and its adjoint 
[Kry05]. As an immediate consequence, the method is symmetric and second-order accurate. 
However, the numerical 'stress test' shows that this integrator is not conjugate symplectic. 

The stress test described in §3 is carried out on LIEMID[EA]. The parameter values and 
initial conditions provided in §3 are used. The time interval of integration is set to be [0, 15000]. 
The energy error obtained on this time interval with a time-stepsize /j = 0.125 is shown in 
Figure 4(a). The experiment was repeated with a time-stepsize h = 0.25 and results are reported 
in Figure 4(b). A systematic drift for the LIEMID[EA] scheme can be observed in both cases. 
The drift appears linear in the time span T and quadratic in the time-stepsize h. We note that the 
drift of LIEMID[EA] is smaller than the drift exhibited by NMB. The time-precision diagrams, 
shown in Figures 5(a) and 5(b) confirm that LIEMID[EA] is second-order accurate. 
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Figure 4: This figure shows the energy error of the Explicit Lie-Midpoint (LIEMID[EA]) and 
velocity Lie-Verlet (VLV) algorithms for the rigid body in the potential energy landscape de- 
fined by (7) for two different timesteps. LIEMID[EA] exhibits a systematic energy drift. On 
the other hand, the energy error of VLV method remains bounded as predicted by theory. The 
initial conditions and parameters used are provided in the text. 
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Figure 5: This figure shows the global error of the Explicit Lie-Midpoint (LIEMID[EA]), Lie- 
Newmark (NMB) and velocity Lie-Verlet (VLV) algorithms. The global error is evaluated 
in configuration and body angular velocity at a physical time of T = 5 for a variety of time- 
stepsizes. We use as a reference solution an integration of (1) using the Matlab function 
ode45 with low tolerance. Observe that all the integrators are second-order accurate. 
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